Load packages

## add 'developer' to packages to be installed from github
x <- c("remotes", "lubridate", "readxl", "pbapply", "viridis", "ggplot2", "kableExtra", "knitr", "formatR", "MASS", "sp", "GGally", "brms", "lme4", "dplyr", "purrr", "forcats", "tidyr", "modelr", "tidybayes", "cowplot", "ggrepel", "posterior", "ggridges", "maRce10/PhenotypeSpace")


aa <- lapply(x, function(y) {
  
  # get pakage name
  pkg <- strsplit(y, "/")[[1]]
  pkg <- pkg[length(pkg)]
  
  # check if installed, if not then install 
  if (!pkg %in% installed.packages()[,"Package"])  {

      if (grepl("/", y))  remotes::install_github(y, force = TRUE) else
    install.packages(y) 
    }

  # load package
  a <- try(require(pkg, character.only = T), silent = T)

  if (!a) remove.packages(pkg)
  })

Functions and global parameters

opts_knit$set(root.dir = "..")

# set evaluation false
opts_chunk$set(fig.width = 10, fig.height = 6, warning = FALSE, message = FALSE, tidy = TRUE)

read_excel_df <- function(...) data.frame(read_excel(...))


# for reading months in english format
sl <- Sys.setlocale(locale = "en_US.UTF-8")


standard_error <- function(x) sd(x)/sqrt(length(x))

cols <- viridis(10, alpha = 0.7)


# print results brms models
summary_brm_model <- function(x, gsub.pattern = NULL, gsub.replacement = NULL, xlab = "Effect size"){
  
  summ <- summary(x)$fixed
  fit <- x$fit  
  betas <- grep("^b_", names(fit@sim$samples[[1]]), value = TRUE)  
  hdis <- t(sapply(betas, function(y)   hdi(fit@sim$samples[[1]][[y]]))
)
  summ_table <- data.frame(summ, hdis, iterations = attr(fit, "stan_args")[[1]]$iter, chains = length(attr(fit, "stan_args")))
  summ_table <- summ_table[rownames(summ_table) != "Intercept", c("Estimate", "Rhat", "Bulk_ESS", "l.95..CI", "u.95..CI", "iterations", "chains")]
  
  summ_table <- as.data.frame(summ_table)  
  summ_table$Rhat <- round(summ_table$Rhat, digits = 3)  
  summ_table$CI_low <- round(unlist(summ_table$l.95..CI), digits = 3)  
  summ_table$CI_high <- round(unlist(summ_table$u.95..CI), digits = 3)  
  summ_table$l.95..CI <- summ_table$u.95..CI <- NULL
  
    out <- lapply(betas, function(y)  data.frame(variable = y, value = sort(fit@sim$samples[[1]][[y]], decreasing = FALSE)))
  
  posteriors <- do.call(rbind, out)
  posteriors <- posteriors[posteriors$variable != "b_Intercept", ]
  
  if (!is.null(gsub.pattern) & !is.null(gsub.replacement))
    posteriors$variable <- gsub(pattern = gsub.pattern, replacement = gsub.replacement, posteriors$variable)
  
   gg <- ggplot(data = posteriors, aes(y = variable, x = value, fill = stat(quantile))) + 
  geom_density_ridges_gradient(quantile_lines = TRUE, quantile_fun = HDInterval::hdi,  vline_linetype = 2) + 
  theme_classic() + 
      labs(x = xlab, y = "Predictor") + 
     xlim(range(c(min(summ_table[ , "CI_low"]), max(summ_table[ , "CI_high"])), 0)) +
  geom_vline(xintercept = 0, col = "red") + 
  scale_fill_manual(values = c("transparent", "lightblue", "transparent"), guide = "none") + ggtitle(x$formula)

  if (!is.null(gsub.pattern) & !is.null(gsub.replacement))
    rownames(summ_table) <- gsub(pattern = gsub.pattern, replacement = gsub.replacement, rownames(summ_table))
   
  summ_table$Rhat <- ifelse(summ_table$Rhat > 1.05, cell_spec(summ_table$Rhat, "html", color ="white", background = "red", bold = TRUE,  font_size = 12),  cell_spec(summ_table$Rhat, "html"))
  
  signif <- summ_table[,"CI_low"] * summ_table[,"CI_high"] > 0
  
  df1 <- kbl(summ_table, row.names = TRUE, escape = FALSE, format = "html", digits = 3)
    
  df1 <- row_spec(kable_input = df1,row =  which(summ_table$CI_low * summ_table$CI_high > 0), background = adjustcolor(cols[9], alpha.f = 0.3))

  df1 <- kable_styling(df1, bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE, font_size = 12)
  
  print(df1)
  
  print(gg)

  }
# read ext sel tab calls
sels <- read.csv("./data/processed/tailored_budgie_calls_sel_tab.csv")

# keep only spectrographic parameters
sels <- sels[, c("sound.files", "selec", "duration", "meanfreq", "sd", "freq.median",
    "freq.IQR", "time.IQR", "skew", "kurt", "sp.ent", "time.ent", "entropy", "meandom",
    "mindom", "maxdom", "dfrange", "modindx", "meanpeakf")]

sels$ID <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][1])
sels$month <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][2])
sels$day <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][3])
sels$year <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][4])

sels$date <- paste(sels$day, substr(sels$month, 0, 3), sels$year, sep = "-")

sels$date <- as.Date(sels$date, format = "%d-%b-%Y")

# acoustic measurements
areas_by_week <- readRDS("./data/processed/acoustic_space_area_by_individual_and_week.RDS")

indiv_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_first_week_by_individual.RDS")

indiv_ovlp$treatment <- factor(indiv_ovlp$treatment, levels = c("Control", "Medium Stress",
    "High Stress"))

group_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_group_by_week.RDS")

# group_ovlp <- do.call(rbind, lapply(unique(group_ovlp$ID), function(x) { X <-
# group_ovlp[group_ovlp$ID == x, ] X$overlap.to.group <- X$overlap.to.group -
# X$overlap.to.group[X$week == min(as.numeric(as.character(X$week)))]
# X$distance.to.group <- X$distance.to.group - X$distance.to.group[X$week ==
# min(as.numeric(as.character(X$week)))] return(X) }))

group_ovlp$treatment <- factor(group_ovlp$treatment, levels = c("Control", "Medium Stress",
    "High Stress"))

Counts per individual

df <- as.data.frame(table(sels$ID))

names(df) <- c("ID", "Sample_size")

df <- df[order(df$Sample_size, decreasing = FALSE), ]

kb <- kable(df, row.names = FALSE)

kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))

print(kb)
ID Sample_size
132YGMM 6
125YGMM 12
160YGHM 16
310YGHM 21
393YGHM 25
303YGHM 37
398WBHM 41
365YGHM 44
399YGLM 46
300YGHM 47
270YGHM 64
407YGHM 97
386WBHM 100
377WWLM 108
367WWNM 119
371YYLM 149
397YGHM 175
378YYLM 195
404WBHM 196
258YGHM 223
389WWLM 230
262YGHM 306
400YGHM 360
388YGLM 377
327YYLM 404
396YBHM 455
403WBHM 566
356WBHM 761
361YGHM 767
402YGHM 849
395WBHM 854
360YGHM 900
390WBHM 935
405YBHM 975
385YBHM 1256
394YBHM 1693

Add metadata

metadat <- read_excel_df("./data/raw/INBREStress_MasterDataSheet_14Nov19.xlsx")

# head(metadat)

sels$ID[sels$ID == "125YGMM"] <- "125YGHM"
sels$ID[sels$ID == "394YBHM"] <- "394WBHM"

# setdiff(sels$ID, metadat$Bird.ID) setdiff(metadat$Bird.ID, sels$ID)

sels$treatment <- sapply(1:nrow(sels), function(x) {

    metadat$Treatment[metadat$Bird.ID == sels$ID[x]][1]

})


sels$treatment.room <- sapply(1:nrow(sels), function(x) {

    metadat$Treatment.Room[metadat$Bird.ID == sels$ID[x]][1]

})


sels$round <- sapply(1:nrow(sels), function(x) {

    metadat$Round[metadat$Bird.ID == sels$ID[x]][1]

})

sels$source.room <- sapply(1:nrow(sels), function(x) {

    metadat$Source.Room[metadat$Bird.ID == sels$ID[x]][1]

})

sels$record.group <- sapply(1:nrow(sels), function(x) {

    metadat$Record.Group[metadat$Bird.ID == sels$ID[x]][1]

})


# add week
out <- lapply(unique(sels$round), function(x) {

    Y <- sels[sels$round == x, ]

    min_date <- min(Y$date)
    week_limits <- min_date + seq(0, 100, by = 7)

    Y$week <- NA
    for (i in 2:length(week_limits)) Y$week[Y$date >= week_limits[i - 1] & Y$date <
        week_limits[i]] <- i - 1

    return(Y)
})

sels <- do.call(rbind, out)


sels$cort.baseline <- sapply(1:nrow(sels), function(x) {

    if (sels$week[x] == 1)
        out <- metadat$D3.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 2)
        out <- metadat$D7.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 3)
        out <- metadat$D14.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 4)
        out <- metadat$D21.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 5)
        out <- metadat$D28.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]

    return(out)
})


sels$cort.stress <- sapply(1:nrow(sels), function(x) {

    if (sels$week[x] == 1)
        out <- metadat$D3.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 2)
        out <- metadat$D7.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 3)
        out <- metadat$D14.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 4)
        out <- metadat$D21.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 5)
        out <- metadat$D28.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]

    return(out)
})

sels$stress.response <- sels$cort.stress - sels$cort.baseline

sels$weight <- sapply(1:nrow(sels), function(x) {

    if (sels$week[x] == 1)
        out <- metadat$D3.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 2)
        out <- metadat$D7.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 3)
        out <- metadat$D14.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 4)
        out <- metadat$D21.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 5)
        out <- metadat$D28.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    return(out)
})


sels$breath.count <- sapply(1:nrow(sels), function(x) {

    if (sels$week[x] == 1)
        out <- metadat$D3.Breath.Count[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 2)
        out <- metadat$D7.Breath.Count[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 3)
        out <- metadat$D14.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 4)
        out <- metadat$D21.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 5)
        out <- metadat$D28.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    return(out)
})


# remove week 5
sels <- sels[sels$week != 5, ]
agg_dat <- aggregate(selec ~ ID + week, data = sels, length)

# compare to week 1 agg_dat$call.count <- sapply(1:nrow(agg_dat), function(x) {
# baseline <- agg_dat$selec[agg_dat$week == 1 & agg_dat$ID == agg_dat$ID[x]] if
# (length(baseline) > 0) change <- agg_dat$selec[x] - baseline else change <-
# agg_dat$selec[x] return(change) } )

# without comparing to week 1
agg_dat$call.count <- sapply(1:nrow(agg_dat), function(x) agg_dat$selec[x])

agg_dat$selec <- NULL

agg_dat$baseline.CORT <- sapply(1:nrow(agg_dat), function(x) {

    baseline <- sels$cort.baseline[sels$week == 1 & sels$ID == agg_dat$ID[x]]
    current <- sels$cort.baseline[sels$week == agg_dat$week[x] & sels$ID == agg_dat$ID[x]]

    if (length(baseline) > 0 & length(current) > 0)
        change <- mean(current) - mean(baseline) else change <- NA

    return(change)
})

agg_dat$stress.response <- sapply(1:nrow(agg_dat), function(x) {

    baseline <- sels$stress.response[sels$week == 1 & sels$ID == agg_dat$ID[x]]
    current <- sels$stress.response[sels$week == agg_dat$week[x] & sels$ID == agg_dat$ID[x]]

    if (length(baseline) > 0 & length(current) > 0)
        change <- mean(current) - mean(baseline) else change <- NA

    return(change)
})


agg_dat$stress.CORT <- sapply(1:nrow(agg_dat), function(x) {

    baseline <- sels$cort.stress[sels$week == 1 & sels$ID == agg_dat$ID[x]]
    current <- sels$cort.stress[sels$week == agg_dat$week[x] & sels$ID == agg_dat$ID[x]]

    if (length(baseline) > 0 & length(current) > 0)
        change <- mean(current) - mean(baseline) else change <- NA

    return(change)
})

agg_dat$weight <- sapply(1:nrow(agg_dat), function(x) {

    baseline <- sels$weight[sels$week == 1 & sels$ID == agg_dat$ID[x]]
    current <- sels$weight[sels$week == agg_dat$week[x] & sels$ID == agg_dat$ID[x]]

    if (length(baseline) > 0 & length(current) > 0)
        change <- mean(current) - mean(baseline) else change <- NA

    return(change)
})

agg_dat$breath.rate <- sapply(1:nrow(agg_dat), function(x) {

    baseline <- sels$breath.count[sels$week == 1 & sels$ID == agg_dat$ID[x]]
    current <- sels$breath.count[sels$week == agg_dat$week[x] & sels$ID == agg_dat$ID[x]]

    if (length(baseline) > 0 & length(current) > 0)
        change <- mean(current) - mean(baseline) else change <- NA

    return(change)
})

agg_dat$acoustic.area <- sapply(1:nrow(agg_dat), function(x) {

    area <- areas_by_week$raref.area[areas_by_week$ID == agg_dat$ID[x] & areas_by_week$week ==
        agg_dat$week[x]]

    if (length(area) < 1)
        area <- NA

    return(area)
})

agg_dat$acoustic.distance <- sapply(1:nrow(agg_dat), function(x) {

    distance <- indiv_ovlp$distance.to.first.week[indiv_ovlp$ID == agg_dat$ID[x] &
        indiv_ovlp$week == agg_dat$week[x]]

    if (length(distance) < 1)
        distance <- NA

    return(distance)
})

agg_dat$acoustic.overlap <- sapply(1:nrow(agg_dat), function(x) {

    overlap <- indiv_ovlp$overlap.to.first.week[indiv_ovlp$ID == agg_dat$ID[x] &
        indiv_ovlp$week == agg_dat$week[x]]

    if (length(overlap) < 1)
        overlap <- NA

    return(overlap)
})

agg_dat$acoustic.overlap.to.group <- sapply(1:nrow(agg_dat), function(x) {

    overlap <- group_ovlp$overlap.to.group[group_ovlp$ID == agg_dat$ID[x] & group_ovlp$week ==
        agg_dat$week[x]]

    if (length(overlap) < 1)
        overlap <- NA

    return(overlap)
})


agg_dat$treatment <- sapply(1:nrow(agg_dat), function(x) unique(sels$treatment[sels$ID ==
    agg_dat$ID[x]]))

agg_dat$round <- sapply(1:nrow(agg_dat), function(x) unique(sels$round[sels$ID ==
    agg_dat$ID[x]]))

Physiological parameters

breath.count <- stack(metadat[, c("D3.Breath.Count", "D7.Breath.Count", "D14.Breath.Count",
    "D21.Breath.Count", "D28.Breath.Count")])

weight <- stack(metadat[, c("D3.Bird.Weight..g.", "D7.Bird.Weight..g.", "D14.Bird.Weight..g.",
    "D21.Bird.Weight..g.", "D28.Bird.Weight..g.")])

cort.stress <- stack(metadat[, c("D3.CORT.Stress", "D7.CORT.Stress", "D14.CORT.Stress",
    "D21.CORT.Stress", "D28.CORT.Stress")])

cort.baseline <- stack(metadat[, c("D3.CORT.Baseline", "D7.CORT.Baseline", "D14.CORT.Baseline",
    "D21.CORT.Baseline", "D28.CORT.Baseline")])


stress <- data.frame(metadat[, c("Bird.ID", "Treatment", "Round", "Treatment.Room")],
    week = breath.count$ind, breath.count = breath.count$values, weight = weight$values,
    cort.stress = cort.stress$values, cort.baseline = cort.baseline$values, stress.response = cort.stress$values -
        cort.baseline$values)

# head(stress)

stress$week <- factor(sapply(strsplit(as.character(stress$week), "\\."), "[[", 1),
    levels = c("D3", "D7", "D14", "D21", "D28"))

stress$day <- as.numeric(gsub("D", "", as.character(stress$week)))
stress$round <- paste("Round", stress$Round)

stress$treatment <- factor(stress$Treatment, levels = c("Control", "Medium Stress",
    "High Stress"))

# remove 5th week
stress <- stress[stress$week != "D28", ]


stress_l <- lapply(stress$Bird.ID, function(x) {
    X <- stress[stress$Bird.ID == x, ]

    X$breath.count <- X$breath.count - X$breath.count[X$week == "D3"]
    X$weight <- X$weight - X$weight[X$week == "D3"]
    X$cort.stress <- X$cort.stress - X$cort.stress[X$week == "D3"]
    X$cort.baseline <- X$cort.baseline - X$cort.baseline[X$week == "D3"]
    X$stress.response <- X$stress.response - X$stress.response[X$week == "D3"]

    return(X)
})

stress <- do.call(rbind, stress_l)

agg_stress <- aggregate(cbind(breath.count, weight, stress.response, cort.baseline) ~
    week + treatment, stress, mean)
agg_stress_se <- aggregate(cbind(breath.count, weight, stress.response, cort.baseline) ~
    week + treatment, stress, standard_error)

names(agg_stress_se) <- paste(names(agg_stress_se), ".se", sep = "")

agg_stress <- cbind(agg_stress, agg_stress_se[, 3:6])

agg_stress$Week <- 1:4

bs <- 22

gg_breath.count <- ggplot(data = agg_stress, aes(x = Week, y = breath.count, fill = treatment)) +
    geom_bar(stat = "identity") + geom_errorbar(aes(ymin = breath.count - breath.count.se,
    ymax = breath.count + breath.count.se), width = 0.1) + scale_fill_viridis_d(begin = 0.1,
    end = 0.9) + facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Mean change in breath\nrate (breaths/min)",
    x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")

gg_weight <- ggplot(data = agg_stress, aes(x = Week, y = weight, fill = treatment)) +
    geom_bar(stat = "identity") + geom_errorbar(aes(ymin = weight - weight.se, ymax = weight +
    weight.se), width = 0.1) + scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment,
    ncol = 3, scale = "fixed") + labs(y = "Mean change in \nweight (grams)", x = "Week") +
    theme_classic(base_size = bs) + theme(legend.position = "none")


gg_cort.baseline <- ggplot(data = agg_stress, aes(x = Week, y = cort.baseline, fill = treatment)) +
    geom_bar(stat = "identity") + geom_errorbar(aes(ymin = cort.baseline - cort.baseline.se,
    ymax = cort.baseline + cort.baseline.se), width = 0.1) + scale_fill_viridis_d(begin = 0.1,
    end = 0.9) + facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Mean change in\nbaseline CORT (ng/mL)",
    x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")


gg_stress.response <- ggplot(data = agg_stress, aes(x = Week, y = stress.response,
    fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = stress.response -
    stress.response.se, ymax = stress.response + stress.response.se), width = 0.1) +
    scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment, ncol = 3,
    scale = "fixed") + labs(y = "Mean change in stress\nresponse CORT (ng/mL)", x = "Week") +
    theme_classic(base_size = bs) + theme(legend.position = "none")

cowplot::plot_grid(gg_breath.count, gg_weight, gg_cort.baseline, gg_stress.response)

Stats

Models: Predicted physio measure ~ treatment + week (continuous) + IndRandom

Variables (Difference from Week 1): weight, BR, baseline CORT, Stress CORT, Stress Response

responses <- c("baseline.CORT", "stress.response", "stress.CORT", "weight", "breath.rate")

predictors <- c("~ treatment + week + (1|ID) + (1|round)")

formulas <- expand.grid(responses = responses, predictors = predictors, stringsAsFactors = FALSE)

vars_to_scale <- c(responses, "week")

# remove week 1
sub_agg_dat <- agg_dat[agg_dat$week != 1, ]

for (i in vars_to_scale) sub_agg_dat[, vars_to_scale] <- scale(sub_agg_dat[, vars_to_scale])

physio_models <- lapply(1:nrow(formulas), function(x) {

    sub_dat <- sub_agg_dat[!is.na(sub_agg_dat[names(sub_agg_dat) == formulas$responses[x]]),
        ]
    sub_dat

    mod <- try(brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
        iter = 20000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9),
        chains = 4), silent = TRUE)
    return(mod)
})


names(physio_models) <- paste(formulas$responses, formulas$predictors)

saveRDS(physio_models, "./data/processed/physiological_response_models.RDS")
physio_models <- readRDS("./data/processed/physiological_response_models.RDS")

ggl <- list()
for (x in 1:length(physio_models)) {
    cat(paste("-  ", names(physio_models)[x], "\n"))
    ggl[[length(ggl) + 1]] <- summary_brm_model(physio_models[[x]], gsub.pattern = "b_treatment|b_",
        gsub.replacement = "")
}
  • baseline.CORT ~ treatment + week + (1|ID) + (1|round)
    Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
    treatmentHighStress 0.763 1.001 10563.66 20000 4 0.114 1.417
    treatmentMediumStress 0.180 1 11354.05 20000 4 -0.511 0.868
    week 0.000 1 30967.77 20000 4 -0.142 0.142
    - stress.response ~ treatment + week + (1|ID) + (1|round)
    Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
    treatmentHighStress -0.493 1.001 8754.839 20000 4 -1.182 0.191
    treatmentMediumStress 0.002 1.001 7335.265 20000 4 -0.768 0.760
    week -0.123 1.001 38086.430 20000 4 -0.255 0.007
    - stress.CORT ~ treatment + week + (1|ID) + (1|round)
    Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
    treatmentHighStress -0.217 1 12610.60 20000 4 -0.908 0.469
    treatmentMediumStress 0.076 1 13662.12 20000 4 -0.697 0.837
    week -0.153 1 37967.80 20000 4 -0.285 -0.020
    - weight ~ treatment + week + (1|ID) + (1|round)
    Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
    treatmentHighStress -0.315 1 10029.64 20000 4 -0.966 0.343
    treatmentMediumStress -0.126 1 10793.59 20000 4 -0.835 0.587
    week -0.344 1 21762.91 20000 4 -0.478 -0.210
    - breath.rate ~ treatment + week + (1|ID) + (1|round)
    Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
    treatmentHighStress 0.145 1.001 21047.468 20000 4 -0.158 0.438
    treatmentMediumStress 0.058 1.001 8702.607 20000 4 -0.275 0.380
    week -0.785 1 35693.191 20000 4 -0.898 -0.670
names(ggl) <- names(physio_models)

 

Barplot and effect sizes side-by-side

cowplot::plot_grid(gg_breath.count + theme_classic(base_size = 10) + theme(legend.position = "none"),
    ggl[[grep("breath", names(ggl))]])

cowplot::plot_grid(gg_weight + theme_classic(base_size = 10) + theme(legend.position = "none"),
    ggl[[grep("weight", names(ggl))]])

cowplot::plot_grid(gg_cort.baseline + theme_classic(base_size = 10) + theme(legend.position = "none"),
    ggl[[grep("base", names(ggl))]])

cowplot::plot_grid(gg_stress.response + theme_classic(base_size = 10) + theme(legend.position = "none"),
    ggl[[grep("response", names(ggl))]])

Takeaways

  • Breath rate decreases gradually with time across after the first week

  • Stress response is higher in “high stress” birds compared to first week

 

Acoustic space projection

t-SNE

scale_param <- scale(sels[, c("duration", "meanfreq", "sd", "freq.median", "freq.IQR",
    "time.IQR", "skew", "kurt", "sp.ent", "time.ent", "entropy", "meandom", "mindom",
    "maxdom", "dfrange", "modindx", "meanpeakf")])

tsne <- Rtsne(scale_param, dims = 2, perplexity = 30, verbose = FALSE, max_iter = 5000)

saveRDS(tsne, "./data/processed/tsne_on_acoustic_parameters_jun_2021.RDS")
tsne <- readRDS("./data/processed/tsne_on_acoustic_parameters_jun_2021.RDS")

Y <- as.data.frame(tsne$Y)
names(Y) <- c("TSNE1", "TSNE2")

sels <- data.frame(sels, Y)

sels$treatment <- factor(sels$treatment, levels = c("Control", "Medium Stress", "High Stress"))

ggplot(sels, aes(x = TSNE1, y = TSNE2, col = as.factor(treatment))) + geom_point() +
    labs(color = "Treatment") + scale_color_viridis_d(alpha = 0.4) + theme_classic(base_size = 25) +
    guides(colour = guide_legend(override.aes = list(size = 10)))

Behavioral parameters

agg_call.count <- aggregate(cbind(call.count, acoustic.overlap.to.group) ~ week +
    treatment, agg_dat, mean)

agg_behav <- aggregate(cbind(acoustic.area, acoustic.overlap) ~ week + treatment,
    agg_dat, mean)

agg_call.count_se <- aggregate(cbind(call.count, acoustic.overlap.to.group) ~ week +
    treatment, agg_dat, standard_error)

agg_behav_se <- aggregate(cbind(acoustic.area, acoustic.overlap) ~ week + treatment,
    agg_dat, standard_error)

agg_behav_se <- merge(agg_call.count_se, agg_behav_se, all = TRUE)

names(agg_behav_se) <- paste(names(agg_behav_se), ".se", sep = "")

agg_behav <- merge(agg_call.count, agg_behav, all = TRUE)

agg_behav <- cbind(agg_behav, agg_behav_se[, 3:6])

bs <- 22

agg_behav$treatment <- factor(agg_behav$treatment, levels = c("Control", "Medium Stress",
    "High Stress"))

gg_call.count <- ggplot(data = agg_behav, aes(x = week, y = call.count, fill = treatment)) +
    geom_bar(stat = "identity") + geom_errorbar(aes(ymin = call.count - call.count.se,
    ymax = call.count + call.count.se), width = 0.1) + scale_fill_viridis_d(begin = 0.1,
    end = 0.9) + facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Vocal output (number of calls)",
    x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")

gg_acoustic.area <- ggplot(data = agg_behav, aes(x = week, y = acoustic.area, fill = treatment)) +
    geom_bar(stat = "identity") + geom_errorbar(aes(ymin = acoustic.area - acoustic.area.se,
    ymax = acoustic.area + acoustic.area.se), width = 0.1) + scale_fill_viridis_d(begin = 0.1,
    end = 0.9) + facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Mean change in \n acoustic area",
    x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")

gg_acoustic.overlap <- ggplot(data = agg_behav, aes(x = week, y = acoustic.overlap,
    fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = acoustic.overlap -
    acoustic.overlap.se, ymax = acoustic.overlap + acoustic.overlap.se), width = 0.1) +
    scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment, ncol = 3,
    scale = "fixed") + labs(y = "Acoustic overlap to self in week 1", x = "Week") +
    theme_classic(base_size = bs) + theme(legend.position = "none")

gg_acoustic.overlap.group <- ggplot(data = agg_behav, aes(x = week, y = acoustic.overlap.to.group,
    fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = acoustic.overlap.to.group -
    acoustic.overlap.to.group.se, ymax = acoustic.overlap.to.group + acoustic.overlap.to.group.se),
    width = 0.1) + scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment,
    ncol = 3, scale = "fixed") + labs(y = "Acoustic overlap to group", x = "Week") +
    theme_classic(base_size = bs) + theme(legend.position = "none")

cowplot::plot_grid(gg_call.count, gg_acoustic.area, gg_acoustic.overlap, gg_acoustic.overlap.group)

Stats

Model: Predicted behavior ~ treatment + week (continuous) + IndRandom

Variables: # calls, Distance moved from self in first week, Overlap to original acoustic space, Match to group repertoire, Maybe overall size of acoustic space

Do as comparison to week one using rarefacted calls and kernel density

responses <- c("call.count", "acoustic.area", "acoustic.overlap", "acoustic.overlap.to.group")

predictors <- c("~ treatment + week + (1|ID) + (1|round)")

formulas <- expand.grid(responses = responses, predictors = predictors, stringsAsFactors = FALSE)

vars_to_scale <- c(responses, "week")


for (i in vars_to_scale) agg_dat[, vars_to_scale] <- scale(agg_dat[, vars_to_scale])

behav_models <- lapply(1:nrow(formulas), function(x) {

    sub_dat <- agg_dat[!is.na(agg_dat[names(agg_dat) == formulas$responses[x]]),
        ]

    # remove week 1
    if (!grepl("count|group", formulas$responses[x]))
        sub_dat <- sub_dat[sub_dat$week != 1, ]

    mod <- try(brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
        iter = 50000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9),
        chains = 4), silent = TRUE)

    return(mod)
})


names(behav_models) <- paste(formulas$responses, formulas$predictors)

saveRDS(behav_models, "./data/processed/behavioral_response_models.RDS")
behav_models <- readRDS("./data/processed/behavioral_response_models.RDS")

ggl <- list()
for (x in 1:length(behav_models)) {

    cat(paste("-  ", names(behav_models)[x], "\n&nbsp;\n---\n&nbsp;"))
    ggl[[length(ggl) + 1]] <- summary_brm_model(behav_models[[x]], gsub.pattern = "b_treatment|b_",
        gsub.replacement = "")
}
  • call.count ~ treatment + week + (1|ID) + (1|round)   —  
    Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
    treatmentHighStress 0.790 1.001 11548.09 50000 4 0.213 1.380
    treatmentMediumStress 0.398 1.001 17599.85 50000 4 -0.182 0.985
    week -0.059 1 102691.50 50000 4 -0.203 0.086
    - acoustic.area ~ treatment + week + (1|ID) + (1|round)   —  
    Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
    treatmentHighStress -0.206 1.001 10719.216 50000 4 -1.059 0.624
    treatmentMediumStress -0.708 1 18329.938 50000 4 -1.567 0.143
    week 0.035 1.003 1444.876 50000 4 -0.094 0.152
    - acoustic.overlap ~ treatment + week + (1|ID) + (1|round)   —  
    Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
    treatmentHighStress 0.910 1 25225.85 50000 4 0.065 1.757
    treatmentMediumStress 0.233 1 25160.95 50000 4 -0.699 1.151
    week -0.169 1 76169.60 50000 4 -0.307 -0.032
    - acoustic.overlap.to.group ~ treatment + week + (1|ID) + (1|round)   —  
    Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
    treatmentHighStress 0.351 1.001 28237.431 50000 4 -0.292 1.015
    treatmentMediumStress 0.380 1.002 3209.631 50000 4 -0.280 1.046
    week -0.223 1.001 41111.148 50000 4 -0.385 -0.057
names(ggl) <- names(behav_models)

 

Barplot and effect sizes side-by-side

cowplot::plot_grid(gg_call.count + theme_classic(base_size = 10) + theme(legend.position = "none"),
    ggl[[grep("count", names(ggl))]])

cowplot::plot_grid(gg_acoustic.area + theme_classic(base_size = 10) + theme(legend.position = "none"),
    ggl[[grep("area", names(ggl))]])

cowplot::plot_grid(gg_acoustic.overlap + theme_classic(base_size = 10) + theme(legend.position = "none"),
    ggl[[grep("overlap ~", names(ggl))]])

cowplot::plot_grid(gg_acoustic.overlap.group + theme_classic(base_size = 10) + theme(legend.position = "none"),
    ggl[[grep("group", names(ggl))]])

Takeaways

  • Higher vocal output in “high stress” birds compared to control

  • Higher acoustic overlap to themselves in week 1 for “high stress” birds compared to control

  • Decrease in overlap to themselves in week 1 across time

 


Session information

## R version 4.1.1 (2021-08-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=es_ES.UTF-8   
##  [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] PhenotypeSpace_0.1.0 warbleR_1.1.27       NatureSounds_1.0.4  
##  [4] seewave_2.2.0        tuneR_1.4.0          ggridges_0.5.3      
##  [7] posterior_1.2.1      ggrepel_0.9.1        cowplot_1.1.1       
## [10] tidybayes_3.0.2      modelr_0.1.8         tidyr_1.2.0         
## [13] forcats_0.5.1        purrr_0.3.4          dplyr_1.0.8         
## [16] lme4_1.1-29          Matrix_1.3-4         brms_2.17.0         
## [19] Rcpp_1.0.8.3         GGally_2.1.2         sp_1.4-6            
## [22] MASS_7.3-54          formatR_1.12         knitr_1.38          
## [25] kableExtra_1.3.4     ggplot2_3.3.5        viridis_0.6.2       
## [28] viridisLite_0.4.0    pbapply_1.5-0        readxl_1.4.0        
## [31] lubridate_1.8.0      remotes_2.4.2       
## 
## loaded via a namespace (and not attached):
##   [1] utf8_1.2.2            tidyselect_1.1.2      fftw_1.0-7           
##   [4] htmlwidgets_1.5.4     grid_4.1.1            munsell_0.5.0        
##   [7] codetools_0.2-18      DT_0.22               miniUI_0.1.1.1       
##  [10] withr_2.5.0           Brobdingnag_1.2-7     spatstat.random_2.2-0
##  [13] colorspace_2.0-3      highr_0.9             rstudioapi_0.13      
##  [16] stats4_4.1.1          dtw_1.22-3            tensor_1.5           
##  [19] bayesplot_1.9.0       labeling_0.4.2        rstan_2.21.5         
##  [22] polyclip_1.10-0       farver_2.1.0          bridgesampling_1.1-2 
##  [25] coda_0.19-4           vctrs_0.4.1           generics_0.1.2       
##  [28] xfun_0.30             R6_2.5.1              markdown_1.1         
##  [31] HDInterval_0.2.2      gamm4_0.2-6           projpred_2.0.2       
##  [34] bitops_1.0-7          spatstat.utils_2.3-0  reshape_0.8.9        
##  [37] assertthat_0.2.1      promises_1.2.0.1      scales_1.2.0         
##  [40] gtable_0.3.0          processx_3.5.3        goftest_1.2-3        
##  [43] rlang_1.0.2           systemfonts_1.0.4     splines_4.1.1        
##  [46] spatstat.geom_2.4-0   broom_0.8.0           checkmate_2.0.0      
##  [49] inline_0.3.19         yaml_2.3.5            reshape2_1.4.4       
##  [52] abind_1.4-5           threejs_0.3.3         crosstalk_1.2.0      
##  [55] backports_1.4.1       httpuv_1.6.5          tensorA_0.36.2       
##  [58] tools_4.1.1           ellipsis_0.3.2        spatstat.core_2.4-2  
##  [61] raster_3.5-15         jquerylib_0.1.4       RColorBrewer_1.1-3   
##  [64] proxy_0.4-26          plyr_1.8.7            base64enc_0.1-3      
##  [67] RCurl_1.98-1.6        ps_1.6.0              prettyunits_1.1.1    
##  [70] rpart_4.1-15          deldir_1.0-6          zoo_1.8-10           
##  [73] cluster_2.1.2         magrittr_2.0.3        ggdist_3.1.1         
##  [76] colourpicker_1.1.1    mvtnorm_1.1-3         matrixStats_0.62.0   
##  [79] shinyjs_2.1.0         mime_0.12             evaluate_0.15        
##  [82] arrayhelpers_1.1-0    xtable_1.8-4          shinystan_2.6.0      
##  [85] gridExtra_2.3         rstantools_2.2.0      compiler_4.1.1       
##  [88] tibble_3.1.6          crayon_1.5.1          minqa_1.2.4          
##  [91] StanHeaders_2.21.0-7  htmltools_0.5.2       mgcv_1.8-36          
##  [94] later_1.3.0           RcppParallel_5.1.5    DBI_1.1.1            
##  [97] boot_1.3-28           permute_0.9-7         cli_3.2.0            
## [100] parallel_4.1.1        igraph_1.3.0          pkgconfig_2.0.3      
## [103] signal_0.7-7          terra_1.5-21          spatstat.sparse_2.1-1
## [106] xml2_1.3.3            svUnit_1.0.6          dygraphs_1.1.1.6     
## [109] svglite_2.1.0         bslib_0.3.1           webshot_0.5.3        
## [112] rvest_1.0.2           stringr_1.4.0         distributional_0.3.0 
## [115] callr_3.7.0           digest_0.6.29         vegan_2.6-2          
## [118] spatstat.data_2.2-0   rmarkdown_2.13        cellranger_1.1.0     
## [121] shiny_1.7.1           gtools_3.9.2          rjson_0.2.21         
## [124] nloptr_2.0.0          lifecycle_1.0.1       nlme_3.1-152         
## [127] jsonlite_1.8.0        fansi_1.0.3           pillar_1.7.0         
## [130] lattice_0.20-44       loo_2.5.1             fastmap_1.1.0        
## [133] httr_1.4.2            pkgbuild_1.3.1        glue_1.6.2           
## [136] xts_0.12.1            shinythemes_1.2.0     stringi_1.7.6        
## [139] sass_0.4.1